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Summary 

A finite difference code has been developed for modeling 
inviscid, unsteady supersonic flow by solution of the com- 
pressible Euler equations. The code uses a deforming grid 
technique to capture the motion of the airfoils and can model 
oscillating cascades with any arbitrary interblade phase 
angle. A flat plate cascade is analyzed, and results are 
compared with results from a small-perturbation theory. 
The results show very good agreement for both the unsteady 
pressure distributions and the integrated force predictions. 
The reason for using the numerical Euler code over a 
small-perturbation theory is the ability to model "rear 
airfoils that have thickness and camber. Sample predictions 
are presented for a section of the rotor of a supersonic 
throughflow compressor designed at NASA Lewis Research 
Center. Preliminary results indicate that two-dimensional, 
flat plate analysis predicts conservative flutter boundaries. 


use a finite difference algorithm to predict the flow field and 
account for these effects. 

The present research focuses on numerical solutions of 
two-dimensional supersonic flow through oscillating cas- 
cades. Numerical methods tend to converge faster in this 
flow regime compared with their subsonic counterparts. The 
Euler equations are solved instead of the Navier-Stokes 
equations since the Reynolds numbers for these flows are 
high. Harmonic pitching motions are prescribed for the 
blade sections for both zero and nonzero interblade phase 
angles. The code uses the deforming grid technique intro- 
duced in references 1 and 2 for convenient specification of 
the periodic boundary conditions from blade to blade. 
Several sample predictions are presented for oscillating 
cascades in this flow regime. This analysis can be applied to 
any fan or compressor designs with supersonic relative Mach 
numbers. 


Introduction 

Recently, there has been interest in designing a turbine 
engine that operates with supersonic axial flow through the 
rotor and stator blade rows. Much research is needed to 
ensure stable operation of this type of design. Understanding 
the aeroelastic behavior and identifying the flutter bound- 
aries is critical. Many methods used for flutter analysis 
depend on the successful predictions of blade loading and 
blade motion. Ideally, this is an interactive process where a 
structural analysis determines the blade motion from the 
blade loading and the loading is determined by the flow 
analysis from the blade motion. Many of the existing 
structural and aerodynamic analysis methods are based on 
flat plate theory ana introduce approximations to model 
thickness and camber of the blades. One alternative is to 


Governing Equations 

A major portion of the present code is based on the 
unsteady solver developed by Sankar and Tang (ref.3) for 
flow past isolated airfoils. This code solves the two- 
dimensional, unsteady, Reynolds-averaged, compressible 
Navier-Stokes equations in strong conservation form on a 
body-fitted moving coordinate system using an ADI proce- 
dure. These equations can be written as 


+ + = + ( 1 ) 


where 
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and p is the fluid density, u and v are the Cartesian compo- 
nents of the fluid velocity, and e is the total energy per unit 
volume. The body- fitted coordinate system is 

related to the Cartesian coordinates by using the following 
transformation: 


S-«*,y,0 

r| = r|(x,y,0 (3) 

T = / 

The Jacobian of the transformation is given by 


^ ^x fly flx^y 


1 

x k y n ~~ x * 


(4) 


and the metrics of the transformation are given by 


z x = J y* flx^-^y* 

Sy = “^*T) ^y = JX l (5) 

Standard central differences were used to compute 
at ^ , y ^ , x n , and y n which were then used to calculate the 
metrics. 

The P and <5 terms in equation (1) are the inviscid terms 
in the ^ and redirections, respectively. The viscous terms, R 
and S, are treated explicitly and can be omitted to give 
solutions of the Euler equations. In fact, all solutions 
presented in this paper are inviscid since the Reynolds 
numbers were relatively high. The Beam-Warming block 
ADI algorithm is used to solve the governing equations. A 
Jameson-type artificial dissipation is added in both direc- 
tions of the computational plane to increase solution sta- 
bility. The solution is second-order accurate in space and 
first- or second-order accurate in time. Further information 
about the algorithm for isolated airfoils can be found in 
reference 3. 


Grid 

A unique feature of the present code is the treatment of 
the grid for oscillating cascades. A method for deforming 
the grid was developed in references 1 and 2 for zero and 
nonzero interblade phase angles. The present study uses the 
same technique. 

The code uses a C-grid generated from Sorenson’s (ref.4) 
GRAPE code, which was modified by Chima (ref.5) for 
improved modeling in turbomachinery problems. One 
C-grid is generated for each blade in the cascade and 
undergoes a local deformation associated with the blade 
motion. The outer boundary of the C-grid is defined by the 
user in the GRAPE code. A deforming grid technique is 
used to locate the position of the grid as a function of time. 
The inner boundary moves with the blade, while the outer 
boundary remains fixed in space. The grid lines connecting 
the inner and outer boundaries are allowed to deform. The 
amount of deformation is a function of the distance away 
from the surface of the airfoil. A weighting function uu is 
defined as 


w(5,fl) = 


s(^fl) 

5 ( ^ ’ fl max ) 


where s is the arc length of a grid line from the airfoil sur- 
face (r|= 1) to some grid point along % - constant, and q max 
is the outer boundary grid line. The grid deformation is 
defined as 


A*iy = «Ay(Ax (7 ) 

Ay^Ufl^Ay;,) (7) 


where A x u and Ay \ i are the spatial differences that would 
exist between successive time steps if the entire grid were 
moved as a rigid body. From equations (6) and (7), we see 
that nodes at the inner boundary (s=0) give w i} - 1, which 
means that the airfoil surface follows the rigid body motion 
of the blade. Conversely, the outer boundary nodes give 
w tJ = 0, and the node positions remain fixed at the initial 
specified locations. The interior nodes shear in space relative 
to the initial grid as w iy varies between 0 and 1. The node 
velocities can be easily found by dividing the grid deforma- 
tion by the time step value. 

Multiple blade computations are made possible by stacking 
the C-grids for each blade and passing information between 
the periodic boundaries. Each C-grid is expanded by one 
grid line in the redirection at the outer boundary to provide 
ghost points for the implementation of the periodic boundary 
conditions. This allows the periodic boundary condition to 
be treated implicitly. The zero interblade phase angle is the 
simplest case for grid generation. Perioaic boundary con- 
ditions are applied across the upper and lower boundaries 
and require a grid for only one blade. However, for 
supersonic axial flow, nonzero interblade phase angles 
require computations around two additional blades for exact 
treatment of the periodic boundary conditions. This 
increases the computational time by a factor of three, but 
provides an exact boundary condition. 

Figure 1 shows a deforming grid for a typical nonzero 
interblade phase angle computation. Multibfade solutions 
for oscillating cascades are done by generating grids for each 
blade and assembling them into a cascade. In this sample 
case, the interblade phase angle a is 180 degrees. A 199 by 
22 C-grid is generated for a typical section near the midspan 
of a supersonic throughflow rotor designed at NASA Lewis 
Research Center. The solver automatically generates the 
grids for the adjacent blades before the unsteady solution 
begins by deforming the mean flow grid for one blade 
through one cycle of oscillation and saving the two grids that 
occur at the desired interblade phase angle. The grids are 
then assembled to give the initial multiblade grid to be used 
for the oscillating cascade solution at uot = 0 (fig. 1(a)). As 
time progresses in the unsteady solution, the grid for each 
blade deforms to model the specified interblade phase angle 
(fig. 1(b)). Notice how the outer boundary of the grid around 
each blade remains fixed in space, while the inner boundary 
follows the motion of the airfoil. An actual run from an 
oscillating cascade typically has a small pitching amplitude 
and does not distort the grid as much as shown in figure 1. 
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Boundary Conditions 

The present solution independently solves the flow 
equations for each grid around a blade and uses periodic 
boundary conditions along the upper and lower boundaries 
to model the cascade effects. Ghost points are assigned at 
the first interior grid line (q = q max - 1) and are used by the 
adjacent grid from the next blade. Although it is tempting 
to use the most current flow information as it becomes 
available from the integration scheme, it is important to only 
use flow information from the same time step across the 
periodic boundaries. This eliminates time inaccuracy due 
to the grid stacking direction of the multiple blade solutions. 
The metric data are also forced to be continuous along the 
periodic boundaries. This procedure essentially makes the 
blade-to-blade periodic boundaries invisible to the flow 
solution. 

In supersonic flow, the domain of the flow field can be 
reduced by inspecting the shock structures through the 
cascade. The disturbances generated inside the Mach cone 
of the leading edge do not influence the flow outside the 
cone, which means that only three blades need to be 
considered in the solutions. For example, consider the 
cascade geometry defined in figure 2, and let the flow condi- 
tions be defined such that the bow shock off the leading edge 
intersects the adjacent blade surfaces. The flow above blade 
2 and below blade 3 will have no influence on the surface 
pressures on blade 1. For this reason, the boundary 
conditions on the upper and lower boundaries of the global 
grid (the grid containing all three C-grids) can be arbitrary. 
The present analysis allows the upper and lower boundaries 
of the global grid to be periodic, which means that only the 
surface pressures from the reference blade will be valid for 
the specified interblade phase angle. 

The inlet conditions are assumed to be uniform by 
specifying the flow density, velocity, flow angle, and pressure. 
The exit flow variables are extrapolated from the interior by 
using a simple first-order model. These boundary conditions 
are valid for supersonic flow at the inlet and exit. Solid wall 
boundary conditions are applied along the airfoil surface, 
and the flow variables are averaged across the slit aft of the 
airfoil. 


Results and Discussion 

Sample predictions for an oscillating cascade in supersonic 
flow have been done for both a cascade of flat plates and a 
midspan section of a proposed rotor design for the super- 
sonic throughflow fan designed at NASA Lewis. There are 
no unsteady experimental data available for comparisons 
with the rotor section predictions. At present, any validation 
of this code is limited to comparisons with a small- 
perturbation theory for flat plates. 


Flat Plate Cascade 

Before an unsteady solution can be obtained for an 
oscillating cascade, a good steady-state solution for the mean 
flow conditions should be calculated. As previously 
described, a C-grid is generated for one blade by using a 
modified GRAPE code. For the present solutions, a 199 by 
22 grid in the and redirections, respectively, is used around 
aflat plate with the thickness-to-chord ratio, x =0.005. The 


leading and trailing edges were rounded to aid in C-grid 
generation, and therefore this is only an approximate 
representation of a flat plate. Three C-grids are then 
assembled to define the geometry for the cascade, as shown 
in figure 3. The present test case was selected from Kielb 
and Ramsey (ref. 6), who simulated the rotor design for the 
supersonic throughflow fan by using Lane’s theory (ref. 7) 
to determine the unsteady pressure distributions. The 
stagger angle y is 28 degrees and the gap-to-chord ratio g/c 
is 0.311. An inviscid, steady-state solution has been gener- 
ated using the present code for M, = 2.61 and P \ = 28 degrees, 
which gives an incidence angle i of 0 degrees. 

Contours of the flow field showing the static pressure are 
shown in figure 4. The shock structure is clearly defined 
showing the reflections of the leading edge Mach waves from 
the surfaces of the adjacent plates. Figure 4(b) shows the 
detail of the pressure contours near the leading edge of the 
plates. There is evidence of a bow shock forward of the 
leading edge, which is expected since the flat plate has finite 
thickness and will demonstrate blunt body flow character- 
istics. This will contribute to deviations from theories that 
use infinitely thin flat plates. Also, there is a slight spreading 
of the shock waves due to the numerical dissipation inherent 
to this type of flow code. It is possible to minimize this 
problem by using a flux vector splitting scheme. In any case, 
the steady-state solution demonstrates the expected flow 
characteristics and should provide a reasonable mean flow 
solution for the oscillating cascade results. The steady-state 
solution required about 1500 iterations with Ax = 0.005 for 
the maximum residual to drop four orders of magnitude, 
which corresponds to about 100 sec of CPU time on the 
CRAY-YMP. This time can be reduced by implementing a 
variable local time step instead of the constant time step used 
in the present solution. Since the focus of the present 
research is the solution of the unsteady Euler equations, little 
time has been spent on accelerating the steady-state solu- 
tions. 

Solutions were done for oscillating flat plates for various 
interblade phase angles and were compared with the 
small- perturbation theory of Lane (ref. 7). The unsteady 
solutions were run with first-order temporal accuracy. The 
surface pressure time histories were recorded and found to 
reach a reasonably periodic solution after three cycles of 
cascade oscillation. A Fourier transform was done on the 
fourth cycle to determine the first harmonic pressure dis- 
tribution relative to the airfoil motion. The pressures were 
normalized by the airfoil pitching amplitude and the phase 
was referenced to the airfoil pitching motion starting at the 
maximum (nose-up) blade angle. Predictions For the 
unsteady pressure difference distributions (A CP) on the flat 
plates are shown in figure 5 for a= 0 and 180 degrees, k= 0.50, 
and a=0.10 degrees. The small amplitude of oscillation is 
used for comparisons with the small-perturbation theory. 
The agreement with Lane’s theory is good for both the real 
and imaginary parts of pressure, except that the locations of 
the shock reflections were predicted about 10 percent 
forward of the Lane’s predictions. The actual locations of 
the shocks are exact for an infinitely thin flat plate, as used 
by Lane’s theory. The finite thickness approximation of the 
flat plate used in the Euler code is expected to be the greatest 
contributor to this error. The wiggles near the leading edge 
of the Euler solution are numerical and are reduced for real 
airfoils. It is possible to adjust the numerical dissipation 
coefficients to minimize these wiggles; however the shocks 
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tend to smear, and the accuracy of the solution becomes 
jeopardized. The force coefficients obtained by integrating 
these pressures are insensitive to these wiggles. 

In an aeroelastic stability analysis, the moment coefficient 
about the elastic axis of tne blade section can be used for 
predicting torsional flutter. If the moment coefficient leads 
the blade motion, the fluid is doing positive work on the blade 
section and flutter is possible if the blade is perturbed by the 
prescribed interblade phase angle and if the structural 
damping cannot control the aerodynamic forces. In the 
oscillating flat plate cascade presented in this paper, an 
unstable moment coefficient is predicted for interblade 

E hase angles near 180 degrees. Figure 6 shows how the time 
istory of the moment coefficient lags the blade motion by 
81.1 degrees (stable) when a=0 and how it leads the plate 
motion by 151.6 degrees (unstable) when a - 180 degrees. 
The moments are calculated about x /c=0.30 (nose-up 
being positive), which also corresponds to the pitching axes. 
Note that the first cycle of the solution contains transients 
from the initialization of the flow field, which was started 
from uniform flow. A Fourier transform was done on the 
fourth cycle to determine the harmonic content of moment 
coefficient. Higher harmonics of the force coefficients were 
found to be small compared to the fundamental frequency. 
Figure 7 summarizes the predictions of the imaginary parts 
of the moment coefficient for various interblade phase 
angles. When the imaginary part of the moment coefiicent 
becomes positive, the moment leads the motion of the flat 
plate and represents an unstable condition. Also shown in 
figure 7 are the predictions using Lane’s theory. The 
agreement is very good over all interblade phase angles in 
spite of the shock locations being different. 

The unsteady solutions use 2.11 x I O ' 5 sec of CPU per 
time step per grid point per blade and required 2047 time 
steps to complete 4.25 cycles of oscillation. In practice, only 
three cycles of oscillation are necessary to reach a periodic 
solution, which corresponds to 400 sec of CRA Y- YMP CPU 
time. 

Rotor Section Cascade 

The reason for using the Euler code over Lane’s theory 
is the ability to model "real" blade sections that have thickness 
and camber. While Lane’s theory offers the most practical 
analysis tool for investigating a wide range of flutter 
parameters, the Euler code can be used to check the effects 
of blade loading. A typical section near the midspan of the 
rotor on a supersonic throughflow fan design being inves- 
tigated at NASA Lewis has been chosen for sample unsteady 
predictions for an oscillating cascade. As in tne flat plate 
results, the steady-state solution is evaluated first to see if 
the mean flow distribution is reasonable. The parameters 
for the flow and cascade geometry are similar to those used 
for the flat plate analysis, except that the airfoils are now 
loaded to conditions near their design point: g/c= 0.311, 
Y=28 degrees, Af, = 2.61, and Pi = 36.0 degrees. The grid 
resented in figure 1 is used in the present investigation and 
as the same dimensions and parameters described for the 
grid used in the flat plate analysis. The steady-state pre- 
dictions for pressure contours and Mach contours are 
presented in figures 8 and 9, respectively. The pressure 
contours give a better definition of the shock structure 
through the cascade. The bow shock, the expansion fan on 
the upper surface, and the compression shock on the lower 


surface are all evident near the leading edge, along with their 
respective reflections from adjacent blades. Also, a weak 
"fish-tail" shock is formed on the trailing edge. These 
predictions show all of the expected physics in a supersonic 
cascade. The addition of thickness and camber to the 
solution caused the shocks to move forward from their 
locations in the flat plate analysis, as shown in figure 4. 

Two unsteady solutions for the blades oscillating with 
a=0 degrees and 180 degrees, A:=0.50, and a=0.1 degrees 
are presented in figure 10. The discontinuities near 20- and 
70-percent chord are evidence of the shock structures shown 
in the steady solution. The shock movement due to a pitching 
amplitude of 0.1 degrees was expectedly small. A stability 
analysis was done by running the code for various interblade 
phase angles, as shown in figure 11. The moment coefficients 
are calculated about a position near the midchord. These 
results are representative of predictions with the pitching 
axis at any location along the chord line. The analysis 
remains stable for all interblade phase angles. In fact, the 
magnitudes of the imajginary part of the moment coefficients 
drop by about one order of magnitude when compared with 
the results from the flat plate analysis in figure 7. The trends 
from this preliminary analysis suggest that the flat plate 
analysis of Lane (ref. 7) is conservative for the prediction of 
flutter. However, this does not show that the rotor design 
is free of any flutter problems since (1) only one frequency, 
one Mach number, and one elastic axis were investigated 
with the Euler code; (2) three-dimensional effects have been 
neglected; and (3) no structural model has been included. 
The frequency used in this investigation was chosen as a 
worst case based on the analysis performed by Kielb and 
Ramsey (ref. 6). Since the solution in reference 6 used the 
mean flow from a flat plate analysis, the critical frequency 
is expected to change tor the loaded cascade. Therefore, 
more frequencies need to be investigated at various Mach 
numbers to ensure that the design remains stable. 


Conclusions 

A code has been developed that solves the nonlinear flow 
field for oscillating cascades in supersonic flow and can be 
used in flutter analysis. A finite difference code has been 
developed for modeling compressible, inviscid, unsteady 
supersonic flow by solution of the Euler equations. The code 
uses a deforming grid technique to capture the motion of 
the airfoils and can model oscillating cascades with any 
arbitrary interblade phase angle. A flat plate analysis is done 
for comparisons with a small-perturbation theory. The 
results show very good agreement for both the unsteady 
pressure distributions and the integrated force predictions. 
Sample predictions are presented for a section of the rotor 
on a supersonic throughflow compressor designed at NASA 
Lewis Research Center. Preliminary results indicate that 
two-dimensional, flat plate analysis predicts conservative 
flutter boundaries. 

The code is a tool for modeling "real" blade sections, 
which can have significantly different flow characteristics 
than results from methods that are restricted to flat plate 
geometries. Obviously, further research is needed to validate 
the code, such as detailed experimental data for oscillating 
cascades. It is possible to add a simple two-dimensional 
structural model to evaluate the effects of structural prop- 
erties. Also, an extension to the three-dimensional flow 
equations may be needed to better model the actual 
compressor design. 
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Appendix - Symbols 


C sonic velocity 

C M moment coefficient 

p - p j 

C P pressure coefficient,^-^ 

c blade chord length 

g total energy of the fluid per unit volume 
g/c gap-to-chord ratio 

9m { } imaginary part of { } 

i incidence angle 

J Jacobian of transformation 

k reduced frequency based on semichord, 

M Mach number 

Re Reynolds number based on chord 

de { } real part of { } 

s arc length of a grid line in the redirection 
t time normalized by ~ 

u , v Cartesian velocities normalized by C , 

V total velocity 

w weighting function for grid deformation 


x ,y Cartesian coordinates 
a amplitude of pitching 

(3 flow angle 

Y stagger angle 

A CP pressure difference coefficient, CP- CP + 
r| normal direction of transformed 

coordinate system 

% chordwise direction of transformed 

coordinate system 
p fluid density 

a interblade phase angle 

t time variable 

x ' thickness-to-chord ratio 

(jo airfoil oscillation frequency 

Subscripts: 

1 conditions at the inlet 

2 conditions at the exit 

+ upper and lower surfaces on airfoil, respectively 
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(a) Global grid. 

(b) Enlarged view of leading-edge region. 

Figure 4.— Static pressure contours for flat plate cascade. Inlet Mach number, M,, 2.61; stagger angle, 7, 28°; gap-to-chord ratio, g/c, 0.311; inlet flow 
angle, 0,, 28° (incidence angle, i, 0°), thickness-to-chord ratio, 7' , 0.005. Blue indicates lower pressure; red indicates higher pressure. 
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Figure 7.— Imaginary part of moment coefficient about x 7c = 0.30 versus 
interblade phase angle; Oat plate cascade. Inlet Mach number, M,, 2.61; 
stagger angle, 7 , 28; gap-to-chord ratio, g/c, 0.311; inlet flow angle, 
28 "±0.10°; reduced frequency based on semichord, k, 0.5. 
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(a) Global grid. 

(b) Enlarged view of leading-edge region. 

Static pressure contours for rotor cascade. Inlet Mach number, M,, 2.61; stagger angle, 7 , 28°; gap-to-chord ratio, g/c, 0.311; inlet flow angle. 
Blue indicates lower pressure; red indicates higher pressure. 
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Figure 9.— Mach number contours for rotor cascade. Inlet Mach number, M h 2.61; stagger angle, 7 , 28°; gap-to-chord ratio, g/c, 0.311; inlet flow angle, 
&\ , 36°. Yellow indicates lower Mach number; pink indicates higher Mach number. 


13 





NASA 

National Aeronautics and 
Space Administration 


1 . Report No. 


NASA TM- 102053 
AIAA-89-2805 


4. Title and Subtitle 


Report Documentation Page 


2. Government Accession No. 


Numerical Analysis of Supersonic Flow Through Oscillating Cascade 
Sections by Using a Deforming Grid 


3. Recipient’s Catalog No. 


5. Report Date 


6. Performing Organization Code 


7. Author(s) 

Dennis L. Huff and T.S.R. Reddy 


9. Performing Organization Name and Address 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 


8. Performing Organization Report No. 

E-4805 

10. Work Unit No. 

535-03-01 

11. Contract or Grant No. 


13. Type of Report and Period Covered 
Technical Memorandum 


1 14. Sponsoring Agency Code 


15. Supplementary Notes 

Prepared for the 25th Joint Propulsion Conference cosponsored by the AIAA, ASME, SAE, and ASEE, 
Monterey, California, July 10-12, 1989. Dennis L. Huff, NASA Lewis Research Center; T.S.R. Reddy, The 
University of Toledo, Toledo, Ohio 43606 and NASA Resident Research Associate. 

16. Abstract 

A finite difference code has been developed for modeling inviscid, unsteady supersonic flow by solution of the 
compressible Euler equations. The code uses a deforming grid technique to capture the motion of the airfoils and 
can model oscillating cascades with any arbitrary interblade phase angle. A flat plate cascade is analyzed, and 
results are compared with results from a small-perturbation theory. The results show very good agreement for 
both the unsteady pressure distributions and the integrated force predictions. The reason for using the numerical 
Euler code over a small-perturbation theory is the ability to model “real” airfoils that have thickness and 
camber. Sample predictions are presented for a section of the rotor on a supersonic throughflow compressor 
designed at NASA Lewis Research Center. Preliminary results indicate that two-dimensional, flat plate analysis 
predicts conservative flutter boundaries. 


17. Key Words {Suggested by Author(s)) 

Oscillating cascades; Supersonic flow; Turbomachinery; 
Euler equations; Numerical analysis 


18. Distribution Statement 

Unclassified - Unlimited 
Subject Category 02 


19. Security Classif. (of this report) 
Unclassified 


20. Security Classif. (of this page) 
Unclassified 


21 . No of pages 
18 


22. Price* 


NASA FORM 1626 OCT 86 


*For sale by the National Technical Information Service, Springfield, Virginia 22161 




